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Abstract 

The problem studied in this paper is ultrasound image reconstruction from frequency -domain mea- 
surements of the scattered field from an object with contrast in attenuation and sound speed. The case 
where the object has uniform but unknown contrast in these properties relative to the background is 
considered. Background clutter is taken into account in a physically realistic manner by considering an 
exact scattering model for randomly located small scatterers that vary in sound speed. The resulting 
statistical characteristics of the interference is incorporated into the imaging solution, which includes 
applying a total-variation minimization based approach where the relative effect of perturbation in 
sound speed to attenuation is included as a parameter. Convex optimization methods provide the basis 
for the reconstruction algorithm. Numerical data for inversion examples are generated by solving the 
discretized Lippman-Schwinger equation for the object and speckle-forming scatterers in the background. 
A statistical model based on the Born approximation is used for reconstruction of the object profile. 
Results are presented for a two dimensional problem in terms of classification performance and compared 
to minimum-^ 2 -norm reconstruction. Classification using the proposed method is shown to be robust down 
to a signal-to-clutter ratio of less than 1 dB. 

I. Introduction 

The last two decades have seen significant interest in the development and evaluation of the lesion 
formation process induced by High Intensity Focused Ultrasound (HIFU) for noninvasive cancer treatment 
fT|-|fT0} • HIFU treatment effectively increases the temperature within an intended region inside the tissue, 
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thus forming a lesion that serves the purpose of eliminating cancerous cells non-invasively. Effective 
monitoring of the lesions is an essential requirement for the treatment to be successful. 

The state of the art in imaging for HIFU treatment is achieved by tracking the temperature changes via 
magnetic resonance imaging. A more convenient and less expensive modality for use in guiding HIFU 
treatment is ultrasound. Recent research in this area includes monitoring HIFU treatment via temperature 



tracking using two dimensional ultrasound 1 1 1 1 . While some studies show that ultrasound backscatter may 



reveal information regarding the temperature distribution following HIFU exposure [12|, it is difficult to 
obtain ultrasound images showing the changes induced by the HIFU process in general, as the backscatter 
properties of the tissue do not change unless cavitation occurs, a process that is typically avoided |T3| , 



|14|. Physics-based approaches |15|, [16| have thus been examined to show the feasibility of HIFU 
monitoring by analyzing traditional ultrasound RF data using wave-based imaging techniques to identify 
the regions that have increased sound speed and acoustic attenuation. 

As a result of HIFU application, the sound speed of the affected tissue has a slight contrast (about 1%) 
but the change is reversible and thus even if imaging is possible based on the sound speed changes fT7) , 



1 18 1, the lesion may no longer be visible once the tissue cools. One property whose contrast changes 
significantly with a rise in temperature is the acoustic attenuation (80-700 %), and the alteration is 
irreversible even after the tissue is cooled JT5J. Therefore, an ultrasound imaging modality based on 



attenuation is potentially suitable to this application. 



Acoustic tomographic imaging is a well-developed field where much of the early work [201, 1 21 1 was 
based on the dual approach of sampling fc-space in circular arcs. In ultrasound imaging for biomedical 
applications there is typically a very limited-view of the region of interest which severely degrades the 
quality of images that can be obtained using diffraction tomography plj. Some of the recent work in 
this area has concentrated on using physics-based models to facilitate the solution of the inverse problem 



1 15), (22) , (23) . The use of the Born model in solving the imaging problem studied here is motivated by 
the simplicity it provides in modeling and computation. In the case of HIFU lesions, the use of the Born 
model is justified partially on physical grounds based on the relatively low contrast with a maximum of 
1% in sound speed and 1.7% in the ratio a p /kf, (the ratio of perturbation in attenuation to the background 
wavenumber) fT3) . While the large size of the lesions relative to the insonifying wavelengh does, strictly 
speaking, violate the Born model; the results of using the Born model to invert data consistent with the 



exact scattering physics will be shown in Sec. IV to be sufficiently accurate in the context of a limited 
view inverse problem to provide suitable results. 

A significant challenge in many sensing contexts is the interference from unknown variations in the 
background that can overwhelm the signal from the object of interest. In the case of ultrasound, sub- 
resolution perturbations to the background properties manifest themselves as speckle in B-mode images. 
There has been significant effort in literature to characterize the formation of speckle using statistical 
methods. An accepted model for the signal magnitude due to coherent integration of the scattered field 
from numerous smaller scatterers is given by the Rayleigh distribution |24[ , p5| . Motivated by this work, 
a Rayleigh distributed model is used in this paper for the amplitude of the speckle forming background 
clutter where the location of the inhomogeneities come from a uniformly random distribution in space. 



The use of total- variation regularization, which is well-known in the inverse problem literature [26|, is 



especially successful when the signal of interest has piecewise smooth features. Under a uniform lesion 



model 1 27 1, HIFU lesions are suitable for imaging with total-variation reconstruction because the binary 
(lesion vs background) classification of the region of interest naturally yields a piecewise-smooth structure. 
In the approach taken here, total-variation reconstruction is used as the objective function of a convex 
optimization problem. The exact model employed in this paper is formed by spatially coincident contrast 
profiles in attenuation and sound speed p"5j . This model allows the characterization of the reduction of 
the number of parameters by a factor of two, due to the relationship between the two profiles. 

The reconstruction method is demonstrated using simulation data generated by solving the Lippman- 
Schwinger equation for the object and the background variation in a discrete setting. While linear addition 



of speckle forming signals |28|-[30| are useful for characterizing the distribution of texture in B-mode 
images (i.e., processed, time domain data), here a fundamentally different approach is taken in that the 
clutter is defined as the spatial distribution of the random, constitutive properties of the background 



medium [24], [25]. This necessitates a first principles approach, described in Sec. Ill where the effect 
of clutter is taken into account by solving the Lippman-Schwinger equation (LSE) in a discretized 
setting to generate the frequency-domain data, based on the spatial distribution of inhomogeneities in 



the background, as will be explained in Sec. IV Solving the LSE enables a more realistic simulation of 
the data contributed by the lesion and background clutter compared to data generated using a matched, 
Born model. Thus, the solution of LSE is performed to include the region of interest where the lesion 



resides, in addition to all background scatterers that reside within the elliptical rings in associated with 
the time-of-flight for the corresponding transmitter-receiver pair. Thermal noise is simulated by linear 
addition of white Gaussian noise to the solution of Lippman-Schwinger equation. The reconstruction 
results are quantified using binary classification performance for the object vs background pixels. Lastly, 
the performance of the proposed method is compared with the minimum-£ 2 -norm reconstruction. 

II. Background 

A. Problem statement 

As motivated in Section [TJ in this paper we are interested in imaging an object that has contrast both 
in sound speed and attenuation with respect to the background. Data are obtained by means of a linear 
array of ultrasonic elements operated in multistatic mode. That is, one element is used as a transmitter 
and then all elements (including the transmit element) are used as receivers to collect backscatter data. 
The system then steps through until all elements have been used to transmit. Further details regarding 



the setting for numerical experiments are specified in Section IV 



B. Mathematical background 

The partial differential equation describing the propagation of time-harmonic acoustic waves is the 



Helmholtz equation [31 1. The Helmholtz equation associated with a constant-density medium is employed 



here per previous work 1 15 1: 



V 2 0(r) + k 2 (r)4){r) = -/(r). (1) 

where <p(r) is the total acoustic field, k(r) = co/c(r) is the wavenumber associated with angular frequency 
oj and sound speed c(r), and /(r) is the source function. For the special case that k(r) = k, a uniform 
background medium, the solution g(r, r') to Eq. ([I} for a point source f(r) = 5{r — r') is defined as the 
free space Green's function, and for 2-D space is given by g(r, r') = (j / 4)H^\k\r —r'\) |3l) . To find the 
total acoustic field in terms of the Green's function, expressing the space-dependent squared-wavenumber 
k 2 (r) as a sum of a constant background component k 2 and a scatter component k 2 (r), one arrives at 

V 2 + k 2 } 0(r) = -f(r) - k 2 s {r)<P{r). (2) 
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Consequently, <p(r) satisfies 

4>{r) = I g(r,r')f(r')dr' + / g(r,r')k 2 (r')(j)(r')dr' (3) 



= <h(r) + [g a <f>](r) (4) 

where the first integral term in Q is recognized as the background field <^>&(r) (the field that would 
be present if there were no inhomogeneities), and Q s represents the last term in Eq. ([3]) in an integral 
operator form that can be interpreted as a mapping from the total field to the scattered field due to k 2 {r). 



This results in the Lippman-Schwinger equation |31|: 



[X-g s ]0(r) = ^(r) (5) 

where X is the identity operator, i.e., [X0](r) = <p(r). The scattered field can then be expressed as the 
difference between the total field and background field, so that 

Mr) = M) ~ Mr) (6) 
= {[1-G S }- 1 -l}Mr) (?) 

When the scattering profile is weak, the operator on the right hand side of Eq. Q may be approximated 
by 



M r ) ~ <5sM r ) = J ' g{r,r')k 2 s (r')M r ')dr' 



(8) 



that is known as the (first-order) Born approximation (2TJ, (3T|. 

Let I, m denote the numbers enumerating the elements of an ultrasound array whose elements are used 
to sequentially transmit pulses and simultaneously receive their returns, where I denotes the transmitting 
element and m denotes the receiving element. The measurements in the frequency domain corresponding 
to element pair (Z, m) are characterized by the Born model to relate physical properties of the object to 



scattered field measurements as formulated in 1 15 1, and are given by 

k,m{w) = / k 2 s (r,u)hi(r,u})h m (r,uj)d 3 r + v lm (u) (9) 
Jv 

such that hi(r,ui) and h m (r,uj) are the spatial transfer functions for transmitting element I and receiving 



element m, found by integrating the corresponding Green's function over the element surface; and ui ;Tn (oj) 
is the associated thermal noise term. The scattering component of the wavenumber is given by fl5| 

, 2 , 2uj , . 2u? . , 2lo 2 . , 

M r > w ) = -j — a p {r,u;) o-c p (r) 5-c a (r) (10) 

c b c 6 h eg 

where c p {r) and a p (r,uj) are the perturbation in sound speed and frequency-dependent attenuation profile 
of the object, respectively, and c s (r) is the background variation in sound speed that produces speckle. 
For soft tissue the attenuation increases approximately linearly with frequency [32] and therefore the 
frequency-dependent attenuation term can be expressed as a p (r,Lo) = tf) p (r) ■ uj/(2tt), where ip p (r) is 
expressed in units of Np/(Hz-m) (15). 

We now invoke the assumption that the object of interest, that is the HIFU lesion, is uniform and both 
the sound speed and attenuation have the same spatial distribution. The ratio of the second term on the 



right hand side of Eq. (10 1 to the first term is thus a constant within the object: 



2TTC p (r) 
c 2 bM r )' 



f^=^zrU- (ID 



The scattering contribution can therefore be expressed as: 

ks(r,u) = bP P (r) -jfiip p {r)) r c s (r). (12) 

The equations expressed above in integral forms are typically discretized and solved using digital 
computers, which requires expressing the variables and data in discrete space. In order to achieve a 
faithful characterization of the Born integral in Eq. ((9]), a sufficient spatial sampling rate is essential. 
The integral representing the Born model Eq. ((9]) used for the reconstruction is discretized in space to 
form a matrix- vector representation of the integral. The rows of the matrix A correspond to the pointwise 
multiplication of the spatial transfer functions such that the entry a L>K at row-i and column-K of the 
matrix A is given by 

a L , K = -— — hi L (r K ,uj L )h mL (r K ,L} L ) (13) 
irc b 

where each different (transmitter, receiver, frequency) combination (l t ,m L ,co t ) is assigned a different row 
l = 1, . . . , M; column-K of A corresponds to a different pixel location r K for k = 1, . . . , N covering 



6 



the region of interest where the lesion resides; and A is the grid area. Let \ m in denote the wavelength 
corresponding to the highest frequency at which measurements are obtained. A spatial sampling rate of 
8 samples per wavelength (A m j n ) was employed in order to achieve acceptable characterization of the 
integral. The frequency domain measurement b L due to the lesion, under the Born approximation, can 
now be expressed as 

N 

b t = (l- jfi) a hK ij) p (r K ) + n L (14) 

K=l 

where n t is the total interference originating from two sources, i.e., the scattering from background 
inhomogeneity and the thermal noise v L at the receiver: 

f —2u 2 

n t = — ^-c s (r)hiXr,0J L )h rnL (r- l uj l )d 2 r + v L . (15) 
J c b 

Hence Eq. (|9]) can be approximated in discrete space by a matrix-vector product 

b=(l-jn)Ax + n (16) 

where x = [ipp(ri) . . . V>p(r/v)] T is the vector of contrast values in attenuation depending on the pixel 
location; b = [b\ . . . b]\j] T is the vector of data points; and n = \n\ . . . i7,n] t is the discrete noise vector 
studied in further detail in the next section. The problem to be addressed is the estimation of x and fj, 
from b. 

C. Statistical modeling of speckle and thermal noise 

The two main sources of noise in the data considered here are the sound speed variations in the 
background and the receiver noise. The receiver (thermal) noise is well characterized by white Gaussian 
distribution, however the characterization of the inhomogeneous background is more complex. Here, for 
purposes of inversion, a statistical approach is used for characterizing the inhomogeneous background, 



where a linear contribution of the background sound speed variation is used via the Born model [28], 



1 33 1 for deriving the statistics of the total noise in the measurements while maintaining mathematical 
tractability. 

The two features characterizing the inhomogeneities in the background are the location and the strength 
of each small scatterer. The physical model for each scatterer is based on the accumulation of sub- 
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pixel changes in the background properties. A Rayleigh density is used to model the resulting amplitude 



distribution of such scattering |24|, [ 25 ] . While there is no exact model in literature replicating the contrast 
in sound speed of sub-resolution scatterers known to the authors at the time of writing, the Rayleigh 
distribution was used to model the sound speed increase due to small scatterers in the background. Data 



in Sec. IV are simulated by first generating the location of each background scatterer randomly, and 
then setting the amplitude of the perturbation from a Rayleigh probability density function, which is 
parametrically controlled to yield varying signal-to-clutter ratio (SCR) levels. 

Lastly, again based on the linear model, the total noise contribution from numerous scatterers weighted 
by the corresponding spatial transfer functions of the transmitting and receiving elements displays a normal 
distribution, which is numerically verified in the next subsection. Therefore, a weighted squared-error 
minimization is appropriate due to its maximum likelihood estimator properties in the linear model (34}, 
where the weighting matrix W is given by the inverse square-root of the noise covariance matrix C n . In 
what follows, the first two moments (mean m n and covariance C n ) of the noise are analytically derived 
in terms of the parameters for the medium and the transducer elements, by using the Born approximation, 

— 1/2 

where m n and W = C n are to be used in the reconstruction in the next section. 



In order to study the statistics of the noise n, Eq. ( |15| ) is analyzed for a model with point scatterers 
with a uniform random spatial distribution in the region of interest, the sound speed variation profile can 
be represented by a sum of delta functions 

c s (r) =^ S j<5(r - r { ) (17) 
i=l 

where N s is the number of scatterers, Sj is the random perturbation value with a Rayleigh distribution for 



the i-th scatterer [33 1, and r« is the corresponding random location. By substituting the above expansion 
into Eq. ( p"5j ) and changing the order of summation and integration, the noise can be expressed as 



i = J^hi i:Uli ( r i) h m l ,LuA r i)si + U L (18) 

i=l C b 

= Y J K i {r i )s l + u i (19) 



i=l 



8 



where 

-2uj 2 

c b 

denotes the Born kernel (35j for transmitter-receiver pair l t -m L , frequency Furthermore, by vectorizing 
the variables in Eq. ([19]) for transducer element pairs (l,m), and frequency uj, the notation can be 
simplified as 

/N, \ 



n 



J2K(r i )-s i )+u (21) 



\i=l 



such that n=[n\... n M ] T , v = [v x . . . vm] T , and K(ri) = [K^n) . . . K M {ri)] T ■ 

D. Calculation of noise statistics 

It has been shown that fully developed speckle results in a Rayleigh distributed amplitude (33J. While 
in the model above individual amplitudes of the scatterers in the background are Rayleigh distributed, the 
summation over a large number of sound speed contrast values multiplied by their corresponding weights 
K{ri) result in a jointly normal interference in the measurement vector. This observation is consistent 



with the Rayleigh envelope of the received signal from speckle [24|, [25| due to the relation between 
circular Gaussian and Rayleigh distributions. In order to verify the joint normality of the observation 
vector, 1000 sets of randomly located scatterers were simulated with Rayleigh amplitude for the example 
setting explained in Sec. [IVJ and the joint normality of the observation vector was verified using Royston's 
multivariate normality test [36] implemented in Matlab (37J. The p- value associated with the Royston 
test (the probability of obtaining a test statistic at least as extreme as the one observed, given that the null 
hypothesis is true, i.e., that the generating distribution is normal in this case) was found to be 0.0864. 
For a significance of 0.05, the null hypothesis is not rejected since p > 0.05, thus justifying the use of 
normal distribution in basing the reconstruction method on the use of mean and covariance as described 
in the next section. 

Following the physical model derived above for random scatterers in the background, the mean m n 
and covariance C n of the noise vector n can be found by 



m n = 

i=l 



^E[K( ri ) Sl ] + E[v] (22) 



= N s ■ E[ Si ] ■ E[K( n )} (23) 
= N s -a s J^/2- [ —K( ri )d 2 n (24) 

where we have assumed v is zero mean, a s denotes the area of the region of integration R s . The 
covariance matrix is computed as 

C n = E[nn*} - m n m* n (25) 

N. 

= 2 ElK^SistK^n)} + E[vv*] - m n m n (26) 
i=l 

N s 

= Y,E[s lS *]E [K(n)K*(n)} + C V - m n m* n (27) 

i=l 

where we used the fact that Ks and n belong to independent noise sources and hence are uncorrelated; 
and that the scattering strength Si is Rayleigh distributed where s» and T{ are uncorrelated, with 



E 



2a' s ■ S(i - j) (28) 

where a s is the mode of the Rayleigh probability distribution. 

If the additive Gaussian noise component is taken to be white, so that C u = a^I, then the noise 
covariance matrix is given by: 

C n = 2a 2 s E[K(ri)K*(n)] + C U - m n m* n , (29) 
i=l 



2a 2 s N s [ —Kir^K^r^n + all-mnm* 



(30) 



The integral in Eq. ( |30| ) can then be approximated by summation over a discrete grid on r, G R s . For 
computing the covariance for the examples in the next section, a spatial sampling rate of 8 samples per 
wavelength (X m in) is used in both dimensions. Figure [T] shows the sub-matrix of C n computed for the 



region of interest used in the examples in Sec. IV corresponding with the use of the middle element 



of the array in both transmission and reception, as a function of frequency, normalized by the largest 
diagonal value of the sub-matrix for a v = and a s = 1. 
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Fig. 1: Covariance sub-matrix corresponding to the use of middle element in both transmission and 
reception at different frequencies. Base- 10 logarithm of the absolute value is shown for each entry. 



III. Framework for estimation 

This section provides the technical framework used to estimate the object profile from frequency 
domain measurements of the scattered field. The whitening of the data is directly incorporated into the 

— 1/2 

optimization problem by including the error-weighting (whitening) matrix W = C n for reconstructing 
the original image x from the measurement vector b. Therefore we aim to find the set of parameters x, 
H that best represents the data as specified by the constraints of the problem, expressed by 



minimize 



TV[(l-jfi)x] 
subject to \\W[(l-jii)Ax-b]\\ 2 <e 
x>0 



(31) 



where TV[(l—jfi)x] is the total variation (1-norm of the gradient) [26] of the contrast profile x multiplied 
by (1 — 31*), A* i s tn e scalar defined in Eq. ( fTT| ), and e is the error radius allowed in the reconstruction. 

The optimality of weighted minimum mean-squared error estimation for measurements with additive 
Gaussian noise as the maximum likelihood estimator is well-known J34j. Generalized Tikhonov regular- 
ization is an example where a matrix-weighting (whitening) is applied to the error term [38]. While a 
constrained minimization form [39] is used here, the use of whitening for mean squared error minimization 
is justified through the equivalence of our formulation to that where the error in the predicted data and 



total variation are linearly combined to form the cost function to be minimized |26|. Piecewise-constant 
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functions tend to have lower total variation as the changes in such functions are sparser, and seeking a 



lower total variation thus favors such solutions |26|, [39]. Therefore a total variation based reconstruction 



is appropriate for HIFU lesions [16|, while the contribution of thermal noise and speckle is accounted 
for by the error radius allowed in the measurements. 

The choice of regularization parameters for total variation reconstruction has been the subject of 
literature based on the relation between the total variation and the error in the observations (26j, (39j- 
|42|. According to the discrepancy principle (43j, the error radius is chosen proportional to the expected 



noise level. While the optimal selection of the regularization parameter in the constrained formulation 
is beyond the scope of the current paper, a consistent method is used in setting the parameter e in 
the examples of the next section. For the imaging of HIFU lesions studied here, an estimate for the 
total noise level in the data can be obtained by insonifying the medium prior to the formation of the 
lesion. Specifically, a realization of prior measurements b s from speckle-forming background variation 
and thermal noise is used to obtain e = ||W6 a ||2/2 used in the reconstruction. 

Due to the multiplication of x and fi, ( |3T| ) is not convex in its variables. Nevertheless, for fixed 
values of fj,, the problem is convex in x. Consequently, ([37]) is solved for a fixed \i. The process is 
repeated for a range of discrete values fi = /zi, . . . , \xm and the solution for which the cost function 
</(//, x) = TV[(1 — is minimum is chosen. The performance of the overall estimation is discussed 
in the next section, where the contrast profile is used in a binary classification of the lesion vs background 
pixels. 

IV. Results 

In order to obtain faithful simulation of the physical principles of acoustic scattering, the Lippman- 
Schwinger equation is solved for the full scattering profile including the object and the background 
scatterers that result in speckle. The frequency domain data from the ultrasound elements are simulated 
by solving Eq. Q on a discrete grid with a spatial sampling interval of A m j rt /8. The discrete Lippman- 
Schwinger equation is solved numerically using the Matlab implementation of the GMRES algorithm, 
which iteratively solves a linear equation to yield a minimum residual over the Krylov subspace ]44| . 
The parameters used in calling the function gmres are the maximum number of iterations set to 500 
and the tolerance in approximation set to 10~ 6 . 



12 



-0.5 0.5 1 -1 -0.5 

x (mm) x (mm) 

(a) (b) 



Fig. 2: (a) Real part of the total pressure field computed using the analytic solution given for the object 



interior in the case of an incident plane wave originating from far left [45 ]. (b) Magnitude of the difference 



between the analytical solution [45] and the numerical solution method employed in this paper. 

In order to verify the validity of the solution in discrete space, the pressure field is compared to an 
analytic solution given for the interior of a circular object in the form of a series expansion for the case 



of an incident plane wave [45 1 for the case of an object with 2 mm diameter and 10 m/s contrast in sound 
speed from a background sound speed of 1540 m/s. The result from the discretized Lippman-Schwinger 
equation 4>ls and the analytic solution <po are shown in Fig. [2a] The mean squared error (the relative 
norm of the difference between the two solutions \\4>ls — ^olh/H^olk) i s l ess tnan 0-27 percent. The 



magnitude of the difference between two calculated fields, i.e., \<j>LS ~ 0o|> is shown in Fig. 2b 

Figure [3] shows the setting used to generate data for the reconstruction for a medium with both a target 
and background scatterers. The data are generated by solving Eq. Q, the Lippman-Schwinger equation, 
for the total contrast profile kg(r) that includes the lesion and the background scatterers as defined in 
Eq. ([10]). Data are generated directly in the frequency domain within 2-5 MHz with a step size of 500 
kHz for each transmitter-receiver combination of all 9 elements of a linear array of point sources with 
inter-element distance of 10 mm. Typically there is significant prior information in monitoring the HIFU 
therapy due to the controlled nature of the lesion formation [jT5|, and thus the region of interest can be 
adjusted according to the available prior information regarding the expected location of the lesion. The 
grid is placed at a 4 mm x 4 mm region around the object. In practice, time-domain systems are used 
to collect data and a time-window is applied in order to reduce interference from regions that are not of 
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interest at the data collection stage. In order to incorporate the effect of time-windowing, we include the 
background scatterers that reside within two spatial ellipses defined by the corresponding time-window 
for each transmitting and receiving element pair, based on the time-of-flight and the size of the region of 
interest. A more detailed explanation on the time-space domain relation is provided in the appendix, along 
with the time-window parameters used in the simulation. The perturbation outside the grid is accounted 
for by setting the solution of Eq. (|5]) to include all pixels that fall within the corresponding timing window 
and have a nonzero contrast relative to the background. The object used in this study is an ellipse of 
dimensions 2 x 2.4 mm centered 50 mm away from the middle element of the array. This size is on the 



smaller range of HIFU lesions |46| , | |47| but keeps the problem numerically tractable. The contrast in 
the sound speed inside the ellipse is 10 m/s relative to the background sound speed of 1540 m/s, and the 
contrast in attenuation is 10 Np/(MHz m) relative to the lossless background. The background variation 
is generated in the form of uniformly located point scatterers. 

The number of background scatterers are set equal to the multiplication of the area of the covered 



region by a density of 250 scatterers per cm 2 similar to 1 30 1 (where a volumetric density of 4000 per 
cm 3 is assumed, i.e. 4000 2 / 3 250). The locations of the scatterers are set by generating uniformly 
random coordinates in x and z in the area shown in Fig. [3] Each scatterer of one pixel size (the same 
pixel size used both in the solution of Eq. |5]) on a discrete grid and also the reconstruction) is given 
a random contrast in sound speed with Rayleigh probability density functions with varying modes to 
obtain different levels of signal-to-clutter ratio (SCR). The SCRs are calculated as 101og 10 (||fy|| 2 /||M| 2 ) 
from the vectors b\ and b s that are associated with the measurements due to the presence of lesion-only 
and speckle-only, respectively. For three different modes of Rayleigh amplitude for scatterers, the SCRs 
are 10.8, 5.2, and 0.3 dB. The corresponding values computed for e = ||W6 s ||2/2 from b s , as described 
previously, are 0.113, 0.215, and 0.307 (normalized by the norm ||6||2 of the total signal b in the presence 
of lesion and speckle together), respectively. Following the solution of the Lippman-Schwinger equation 
for the field scattered from the lesion and the background scatterers, thermal noise is added to the resulting 
scattered field measurements where the noise variance is chosen to yield 30 dB SNR in the scattered 
field measurements. 

Figure |4] shows the contrast profile |1 — j\x\x only due to the presence of the lesion. The true value of 
/i is (2-7T • 10)/(1540 2 • 10~ 5 ) « 2.65, while the estimation is performed for a discrete set of candidate fi 
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Fig. 3: Background variation generated as point scatterers within the shown rectangular region, where the 
object is located within the 4x4 mm box centered at z = 0, x = 50 mm, and the transmitting-receiving 
elements are marked with circles on the x = line. 



values on a grid with 0.1 interval between [0.5,10], hence used in the estimation of x per ( f3"Tj ). 



Figures 5a 5b show the reconstruction results for the case with 10.8 dB SCR. In order to quantify the 
performance of the reconstruction, a binary quantization is employed for the classification of the lesion 
and non-lesion pixels in the reconstructed image. For different values of the threshold, the number of 
pixels correctly identified as lesion ucl, and those that are misclassified as lesion uml are counted. 
Both counts are divided by the number of true lesion pixels titlI resulting in an estimated probability of 
detection Pd = I — ncL/nrL for lesion pixels, and a relative false alarm rate rt a = uml/iT'TL, where the 
normalization for rj a is not a probability but expresses the number of the falsely classified background 
pixels relative to the lesion size. For example rj a = 0.05 corresponds to the case where the number of 
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Fig. 4: Contrast of the true object in attenuation in units of Np/(m MHz), multiplied by |1 — 



misclassified background pixels is 1/20 of the total number of pixels in the true lesion. The performance 
curves for the cases with different SCRs are plotted in Fig. [6] obtained by setting the threshold as integer 
multiples of 1/100 of the peak value, shown in the interval 10 -2 < r/ a < 1. For a relative false alarm 
rate of rj a = 0.05, the probability of detection is near 0.98, 0.93, and 0.84 for the cases with 10.8, 
5.2, and 0.3 dB SCR, respectively. The estimates fi were 1.3, 1.2, and 1.1, respectively, compared to the 
actual value of 2.65. 

To validate that estimates of /j, improve in cases with lower interference, the reconstruction is performed 
for the same sized object with lower contrast in sound speed (1 m/s) and attenuation (1 Np/(m MHz)) for 
an SCR of 11.8 dB, where the error in the Born model is smaller due to the lower contrast relative to the 
background. Figures [7] and [8] show the associated reconstruction results and the classification performance. 
In this instance, the estimate of /i is 2.3, significantly closer to the true value 2.65, and the detection 
probability was pd ~ 0.99 for rf a = 0.05. Therefore, while the estimation of \i is affected by the total 
interference in the measurements (model error plus noise), the classification performance curves in Fig. [6] 
indicate that the ability to estimate the boundary of the lesion is robust to estimation errors in fi. 

Lastly, we compare the performance of the proposed method to the case where the £ 2 -norm of the 
variable x is used as the objective function instead of total variation. The use of £ 2 -norm as a penalty 



function added to the residual is well-known as Tikhonov regularization [38], [40 1, which improves the 
conditioning of the problem to enable a numerical solution. In order to maintain a similar setting in 



comparison to the application of our method, an optimization form, similar to ( 3 1 1 is used, except that fi 
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(b) Binarized reconstruction and borders of the actual lesion 



Fig. 5: (a) Contrast of the estimated object in attenuation in units of Np/(m MHz), multiplied by the 



magnitude |1 — jfi\ for each case; (b) binarized (lesion vs non-lesion) reconstruction where the edges of 



the actual lesion are outlined with the solid black line. 
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Fig. 6: Binary classification performance curves (detection probability vs relative false alarm) for SCRs 
of 10.8, 5.2, and 0.3 dB SCR 



is assumed known for the objective function and data constraint. Specifically, the following optimization 
problem is solved 

minimize II x II 2 

X 

subject to ||W[(1 - jfi)Ax - b)\\ 2 < e (32) 
x > 
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(a) Reconstruction of object with 1/10 the contrast for 11.8 dB SCR (b) Binarized reconstruction and borders of the actual lesion 
Fig. 7: 



(a) Contrast of the estimated object in attenuation in units of Np/(m MHz), multiplied by the 



magnitude |1 — jfi\ for the object contrasting 1 m/s and 1 Np/(m MHz) in sound speed and attenuation, 



respectively; (b) binarized (lesion vs non- lesion) reconstruction and the edges of the actual lesion. 




Fig. 8: Binary classification performance curve for lesion with 1/10 of the contrast and 11.8 dB SCR. 



where the true value \i = 2.65 is prescribed, and all other parameters are the same as in ( |31j ). Figures 9a 
shows that, despite the use of the true value of \i, the minimum-£ 2 -norm method is not able to reconstruct 
the object as well as the TV algorithm: the object is not dark inside and it is surrounded by ripples. For 
comparison purposes, Fig. [9b] shows the resulting classification performance from these reconstructions. 
The performance curve is significantly worse with a failure of the classification pd < 0.13 for r/ = 0.05 
compared to pd « 0.98, pd ~ 0.94, and pd ~ 0.81, for SCRs of 10.8, 5.2, and 0.3, respectively, as seen 
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(a) Min-^ 2 -norm reconstruction for 10.8 dB SCR 
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(b) Binary classification performance 



Fig. 9: (a) Contrast of the object in attenuation for minimum-^ -norm estimation, multiplied by the 



magnitude |1 —jfi\ at 10.8 dB SCR; (b) Detection probability vs relative false alarm for the three cases. 



in Fig. [6] for the TV based method. In comparison, the performance of the TV based is even better when 
the true value fj, = 2.65 is prescribed, resulting in p^, 0.99, p,i ~ 0.95, and pa ~ 0.83, for SCRs of 
10.8, 5.2, and 0.3, respectively. 

V. Conclusions and future work 

In this paper a method was presented to image a lesion with uniform sound speed and attenuation 
profile in 2-D. The methods used in this paper can be extended to problems in three-dimensional case 
with sufficient engineering effort for solving the convex optimization problems in higher dimensions. 
The method used for reconstruction is based on the Born model for frequency-domain measurements 
of the backscattered ultrasound signal from a linear array of elements. By relating the scattered field 
measurements to the contrast profiles of the object in sound speed and attenuation, where one profile is 
expressed as a scalar times the other, the dimension of the estimation problem is reduced by half. 

A model for the speckle forming clutter was employed using a random distribution of point sources 
where the scattering amplitudes were chosen so that the statistics matched a Rayleigh distribution. The 
Born approximation was then used to derive the statistics of the total noise in the observations. In order 
to take advantage of the structure of the interference due to the inhomogeneous background consisting 
of randomly located small scatterers, we performed a covariance analysis and applied whitening to the 
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measurement model. 

For the numerical experiments, frequency-domain data were generated by solving the Lippman-Schwinger 
for the lesion, where the background scatterers were also included in the solution of the Lippman- 
Schwinger equation. The validity of the numerical solution to Lippman-Schwinger equation was checked 
by comparison to analytical expressions for a case previously studied in literature. Thermal noise in the 
measurements was simulated by linearly adding zero-mean white Gaussian noise, while speckle generating 
clutter was accounted for directly as perturbation in sound speed in the solution of Lippman-Schwinger 
equation. 

For solving the problems ([31) and (|32"|) we used CVX, a package for specifying and solving convex 



programs B9J, [50 1. For larger problems with higher computational load, it is possible to implement 



problem-specific solutions to the convex optimization problem ([31]) to avoid the overhead of using general 



purpose convex solvers [40|, |51|. 



Our study of the classification performance shows the feasibility of the reconstruction method down to 
a signal-to-clutter ratio of less than 1 dB for identifying the lesion and background pixels. The detection 
rate of the TV based reconstruction method generally exceeded 80%. In comparison, a minimum-£ 2 -norm 
reconstruction under similar constraints even when the true value of the scalar parameter is assumed to 
be known, had a detection rate of less than 15%. These results suggest that TV based reconstruction is 
inherently more suited to the piecewise smooth nature of the HIFU lesions, as well as that a TV approach 
to detecting HIFU lesions can make detection of lesions by ultrasound practical. 

Appendix 
Time windowing of measurements 

As mentioned in the referring text above, time-domain systems are used in practice in order to obtain 
ultrasound measurements. Time windows are typically applied to avoid large signal returns from nearby 
scatterers. In our example application, the expected region for the object (HIFU lesion) location is known 
by design, and therefore application of a time-window is useful for avoiding excess interference that are 
due to background scatterers outside the region of interest. 

While simulating the application of a time-window involves solving the LSE for the entire domain 
of scatterers for all frequencies of interest, finding the corresponding time-domain sequence by inverse 
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Fourier transformation, and then applying a time-window to the resulting sequence; application of a 
spatial window involves a Born approximation for scatterers that are outside a close time -proximity of 
the object of interest by ignoring the secondary interactions between such scatterers. Due to the small 
size and scattering strength of the sub-resolution scatterers used in our numerical experiments, we adopt 
this approach in the simulations to make the simulation feasible. Note that all secondary interactions 
within the spatial region associated with the time-window for a given element pair are accounted for in 
the solution of Eq. ([5]), including those of the lesion pixels and the background scatterers, both between 
themselves and each other. Below arguments are thus applicable under these conditions. 

For a given medium with approximately constant sound speed, the time-of-flight of a pulse is determined 
by the distance it travels in space. For a point element, the scattering location for a pulse emitted at t = 
and returning at t = ti resides on a circle with radius r\ = c ■ t\/2 where c denotes the speed of 
sound. Therefore, in a homogeneous background, a time-window applied between [ii,<2] is associated 
with an annular region whose distance to the point element remains between [r\, r-z\, where r2 = c-t^ll. 
For the examples above, the spatial window at a distance between [47, 53] mm from the center element 
corresponds to a time-window of [61.04,68.83] fis for c = 1540 m/s. For a pulse with 3.5 MHz center 
frequency and 3 MHz bandwidth, the duration is on the order of 0.3 fis and thus a time-window width 
of T = 7.8 [is covers an interval of sufficient width for processing. 

The same window width is used for each transmitter-receiver pair where the center of the time-window 
is selected as the sum of the distances of the center of the region of interest to the transmitting and 
receiving elements divided by c. In spatial terms, we use only those scatterers that reside between the 
two ellipses with the two common focal points (location of transmitting and receiving elements) and radii 
Ti and r a , where r j = r c — Tc/4 and r a = r c + Tc/A where r c is the sum of distances of the center of 
the region of interest to the two focal points. Figure [10] shows the spatial regions corresponding to the 
time-domain windows for element-pairs (l,m) = (5,5) (middle element used in both transmission and 
reception), (2, 2) (second element from left used in both transmission and reception), and (3, 9) (the third 
element from left and the rightmost element used in transmission and reception, respectively), where the 
elements are numbered from 1 to 9 (from left to right). Note that the the regions remain unaltered when 
the transmitting and receiving elements are swapped, i.e., element pair (l,m) and (m, I) use the same 
interval for time-windowing in the receiving element and have identical spatial regions in association. 
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z (mm) 



Fig. 10: Elliptical rings showing the borders of the regions associated with the timing window for the 
corresponding element pairs of (l,m) = (5,5) (solid), (2,2) (dashed), and (3,9) (dotted). The elements 
are marked with circles and have coordinates of x = 0, z = {— 40, — 30, 40} in millimeters. The 
region of interest (shown with a box) is located at the intersection of the annuli. 
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